Visualizing Spatial Data in R

Kristin Broms, Neptune & Co.

October 18, 2017

Outline

As you can see, “spatial data” has come to mean much more than creating a map! Visualizing spatial data is an active area of software development, with many exciting packages and functions becoming available on a regular basis.

Static map- shapefiles

library(rgdal)
streams <- readOGR(dsn = "../data/R_GIS_Layers/", 
                   layer = "Cove_Drainage_WGS84")
## OGR data source with driver: ESRI Shapefile 
## Source: "../data/R_GIS_Layers/", layer: "Cove_Drainage_WGS84"
## with 49 features
## It has 9 fields
summary(streams)
## Object of class SpatialLinesDataFrame
## Coordinates:
##          min        max
## x -109.35383 -109.12624
## y   36.46573   36.94164
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##     OBJECTID          ComID               GNIS_ID       GNIS_Name 
##  Min.   :     0   Min.   :        0   00003435: 2   Cove Wash: 2  
##  1st Qu.: 32992   1st Qu.:103664065   NA's    :47   NA's     :47  
##  Median : 65593   Median :103665843                               
##  Mean   : 68116   Mean   :101439114                               
##  3rd Qu.:105309   3rd Qu.:126660550                               
##  Max.   :136421   Max.   :126669825                               
##                                                                   
##           ReachCode      FType         FCode        Source  
##  14080204001554: 2   Min.   :460   Min.   :    0   DOQQ: 1  
##  14080204001982: 2   1st Qu.:460   1st Qu.:46003   NHDH:48  
##  14080105001404: 1   Median :460   Median :46003            
##  14080105001406: 1   Mean   :462   Mean   :45264            
##  14080105001455: 1   3rd Qu.:460   3rd Qu.:46003            
##  (Other)       :41   Max.   :558   Max.   :55800            
##  NA's          : 1                                          
##                Name   
##  Middle 1        : 2  
##  Cove Wash       : 1  
##  Cove Wash Middle: 1  
##  Cove Wash North : 1  
##  Cove Wash South : 1  
##  (Other)         :17  
##  NA's            :26

Static map- shapefiles

plot(streams, col="blue")

## note that plots using base R may be distorted.

Static map- shapefiles

library(ggplot2)
streams_gg <- fortify(streams)  ## or use broom::tidy()
# str(streams_gg)
ggplot() + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
                color = "blue", size = 1.5)

## note that without a coordinates correction, the plot may be distorted.

Static map- satellite maps

library(ggmap)
devtools::install_github("hadley/ggplot2@v2.2.0")
library(ggplot2)
myLocation <- c(lon = -109.2279, lat = 36.545) 
myMap <- get_map(location = myLocation, source = "google", 
                 maptype = "satellite", crop = FALSE, zoom = 13)
g <- ggmap(myMap, darken = c(0.25, "white"))
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5)
ggsave("shp_and_ggmap.png", dpi = 72) 

Static maps - with data and formatting

library(dplyr)
alum_data <- plot_data %>% filter(Analyte == "ALUMINUM")
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5) +
  geom_point(data=alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             shape = 21, alpha = 0.9) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
  scale_size(range = c(2, 8)) +
  guides(size = FALSE) +
  coord_equal() 
ggsave("static_map.png", dpi = 72) 
head(alum_data)
##    Analyte    Result Units Latitude Longitude
## 1 ALUMINUM  9120.859 mg/kg 36.59625 -109.1736
## 2 ALUMINUM  7429.408 mg/kg 36.59625 -109.1736
## 3 ALUMINUM  1657.721 mg/kg 36.61342 -109.1356
## 4 ALUMINUM  1709.072 mg/kg 36.61342 -109.1356
## 5 ALUMINUM  5867.922 mg/kg 36.56392 -109.2120
## 6 ALUMINUM 46050.939 mg/kg 36.56392 -109.2120

Static map - with data and formatting

Using R in place of ArcGIS

3D maps

Lollipop example, using scatterplot3d

library(scatterplot3d)
with(alum_data, {
  lollipop_plot <- scatterplot3d(streams_gg$long, streams_gg$lat, 
                                 rep(0, length(streams_gg$long)),  # x y and z axis
                              color = "blue", pch = 16, type = "p",
                              angle = 120,            
                              scale.y = 1.75,           
                              main = "Example Lollipop plot",
                              xlab = "Longitude",
                              ylab = "Latitude",
                              zlab = "Concentration values",
                              xlim = c(-109.273, -109.15),
                              ylim = c(36.50, 36.59),
                              zlim = range(0, 46000),
                              grid = TRUE,
                              box = FALSE)
  
  # add the legend
  legend("topright", inset = 0.07,      
         bty = "n", 
         title = "Type of Results",
         c("High","Low"), 
         fill = c("#1B9E77", "#D95F02"))
  
  # add the lollipop points
  lollipop_plot$points3d(Longitude, Latitude, Result,        
                      col = "#282830", 
                      pch = 21, 
                      bg = ifelse(Result > 10000, "#1B9E77", "#D95F02"), 
                      lwd = 1,        
                      type = "h",      
                      cex = (Result / 50000) + 1)
})

Lollipop example

Lollipop example

lollipop plot

lollipop plot

(Generated using a modified version of the scatterplot3d function: https://github.com/USEPA/R-User-Group/tree/master/contributedCode/scatterplot3dMap)

Reference: Beaulieu, et al. 2016. Estimates of reservoir methane emissions based on a spatially balanced probabilistic-survey. Limnology and Oceanography, 61: S27-S40.

3D map example (interactive)

library(rgl)
library(dplyr)
plot_3d <- with(alum_data, 
  plot3d(Longitude, Latitude, log(Result),        # x y and z axis
         col = ifelse(Result > 10000, "#1B9E77", "#D95F02"), size = 5,
         type = "p"))
rglwidget(elementId = "plot3drgl") # to show in presentation

3D map- plotly

## plotly recommends the development version of ggplot2, but will also work with the
## latest version of the package, version 2.2.1 
# devtools::install_github('hadley/ggplot2')
# library(ggplot2)
library(devtools)
library(dplyr)
library(plotly)

p <- plot_ly(alum_data, 
             x = ~Longitude, y = ~Latitude, z = ~Result,
             marker = list(color = ~Result, 
                      colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
  add_markers() 

3D map- plotly

p

Interactive maps- plotly

interactive_map <- g +  
  geom_point(data = alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             color = "#000000", shape = 21, alpha = 0.8) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728",
                      guide = "legend") +
  scale_size(range = c(2, 8)) +
  coord_equal() +
  ggtitle("ALUMINUM")

Interactive maps- plotly

ggplotly(interactive_map + theme_void(), filename="plotly")

Interactive maps cont’d

Animated maps

library(animation)
ani.record(reset = TRUE)
for(a in unique(plot_data$Analyte)) { 
  plot_data_a <- subset(plot_data, Analyte == a)
  animated_maps <- g + 
    geom_path(data = streams_gg, 
              aes(x = long, y = lat, group = group),
              size = 1.5, color = 'lightblue') +
    geom_point(data = plot_data_a, 
               aes(x = Longitude, y=Latitude, size = Result, fill = Result),
               shape = 21, alpha = 0.8) +
    scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
    scale_size(range = c(2, 8)) +
    guides(size = FALSE) +
    coord_equal() +
    ggtitle(a)
  
  print(animated_maps)
  ani.record()
}
oopts = ani.options(interval = 1)
ani.replay()
saveHTML(ani.replay(), img.name = "animation_plot", 
         htmlfile = "animation_5metals.html", ani.width=600, ani.height=400)

Animated maps

Spatial data exploration- Moran’s I

# Show neighbors on map of subbasins
library(sp)
class(lulc.sp1)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(lulc.sp1, zcol = "PCTHERB", col = "black", 
       main = "Percent Herb")

Spatial data exploration- Moran’s I

library(spdep)

## determine who shares a border
ctchmt.nb1 <- poly2nb(lulc.sp1, queen = TRUE)
# put neighbors into list
ctchmt.nbwts.list <- nb2listw(ctchmt.nb1, style = "B")

moran.test(lulc.sp1$PCTHERB, ctchmt.nbwts.list)
## 
##  Moran I test under randomisation
## 
## data:  lulc.sp1$PCTHERB  
## weights: ctchmt.nbwts.list  
## 
## Moran I statistic standard deviate = 2.9826, p-value = 0.001429
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.212182899      -0.016949153       0.005901866

Spatial data exploration- Moran’s I cont’d

par(mar=c(5, 5, 2, 2))
moran.plot(lulc.sp1$PCTHERB, listw = ctchmt.nbwts.list, 
           labels = lulc.df1$PCTHERB,
           xlab = "PCTHERB", ylab = "Spatially Lagged PCTHERB",
           main = "Moran Plot for PCT HERB")

Additional Packages

Appendix: Interactive network graph

library(magrittr)  ## for %<>% function
library(dplyr)  
library(Rgraphviz)  ## to convert data to necessary data types
library(visNetwork)  ## to create the interactive network graph
network_interactive <- function(g, nodes, coords){
  ## first, switch graph construction to data frames:
  nNodes <- length(nodes)
  nodes_df <- matrix(NA, nrow=nNodes, 4)  # save first 4 attributes
  for (i in 1:nNodes){
    nodes_df[i, ] <- unlist(jsonlite::fromJSON(nodes[[i]]))[1:4]
  }
  nodes_df <- as.data.frame(nodes_df, stringsAsFactors=FALSE)
  
  ## because edges uses the names, need id to equal names of nodes, not numbers
  names(nodes_df)[1:3] <- c("old_id", "id", "type")
  
  ## and id column needs to go first: (only for igraph package, not visNetwork)
  ## igraph packge was used to help decide node coordinates
  nodes_df <- nodes_df[, c('id', 'old_id', 'type')]

  g_edges <- Rgraphviz::buildEdgeList(g)
  nEdges <- length(g_edges)
  edges_df <- NULL
  for (i in 1:nEdges){
    tmp <- c(from=g_edges[[i]]@from,
             to=g_edges[[i]]@to,
             unlist(g_edges[[i]]@attrs))
    edges_df %<>% bind_rows(tmp)
  }
  edges_df$weight <- as.numeric(edges_df$weight)
  
  ## Add layer info/ color column
  nodes_df$col_layer <- 
    ifelse(nodes_df$id %in%  grep("MainInput", nodes_df$id, value=T), 1, # 
           ifelse(nodes_df$id %in% grep("OtherInput", nodes_df$id, value=T), 3,
           ifelse(nodes_df$id %in%  grep("Midvalue", nodes_df$id, value=T),  4,
           ifelse(nodes_df$id %in%  grep("MidEqn", nodes_df$id, value=T),  5,
                  6))))
  
  ##  add coordinates to determine layout of each node:
  nodes_df <- full_join(nodes_df, coords)

  ## specify colors to use for each col_layer
  graph_colors <- c("gold", "darkorange", "tomato", 
                    "palegreen", "seagreen", "royalblue")
  # frame = borders
  frame_colors <- c("gold3", "darkorange3", "tomato4", 
                    "palegreen3", "seagreen4", "royalblue4")
  ## (visNetwork doesn't like color names with numbers)
  frame_colors_rgb <- rgb(t(col2rgb(frame_colors)), maxColorValue = 255)

  ## Add attributes so that the graph looks good:
  nodes_df$shape <- "ellipse"
  nodes_df$color.background <- graph_colors[nodes_df$col_layer] 
  nodes_df$color.border <- frame_colors_rgb[nodes_df$col_layer]
  nodes_df$color.highlight.border <- "darkred"
  
  edges_df$arrows <- "to"  # draw arrows on the edges.
  edges_df$color.highlight <- "black"
  
  network <- visNetwork::visNetwork(nodes_df, edges_df, width="100%", physics = TRUE) %>%
    visNetwork::visIgraphLayout() %>%
    visNetwork::visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE,
               manipulation = FALSE) %>%
    visNetwork::visEdges(color = "slategray", smooth = list(roundess = 1)) %>%
    visNetwork::visNodes(font = list(size = 22), shape="ellipse")

  return(network)
}

Appendix: Interactive network graph

network_interactive(g = networkDag, nodes = networkNodes, coords = networkCoords)

Appendix: Color

Appendix: Color

par(mar=rep(0, 4))  ## change margins so plot fills entire space
## default colors:
plot(1:8, pch=16, cex=5, col=1:8)

## Defining a new palette
new_palette <- c("darkred", "chartreuse", "turquoise", "purple",  
                 "gray45", "plum", "black", "#F08800")
palette(new_palette)
plot(1:8, pch=16, cex=5, col=1:8)

palette("default")  ## return palette to default colors

Version Control

sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] visNetwork_2.0.1    Rgraphviz_2.18.0    graph_1.52.0       
##  [4] BiocGenerics_0.20.0 magrittr_1.5        spdep_0.6-15       
##  [7] Matrix_1.2-11       bindrcpp_0.2        plotly_4.7.1       
## [10] devtools_1.13.3     dplyr_0.7.3         rgl_0.98.1         
## [13] ggplot2_2.2.1.9000  rgdal_1.2-11        sp_1.2-5           
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.3.1        maps_3.2.0        tidyr_0.7.1      
##  [4] jsonlite_1.5      viridisLite_0.2.0 splines_3.3.3    
##  [7] gtools_3.5.0      shiny_1.0.5       assertthat_0.2.0 
## [10] expm_0.999-2      stats4_3.3.3      yaml_2.1.14      
## [13] LearnBayes_2.15   backports_1.1.0   lattice_0.20-35  
## [16] glue_1.1.1        digest_0.6.12     colorspace_1.3-2 
## [19] htmltools_0.3.6   httpuv_1.3.5      plyr_1.8.4       
## [22] pkgconfig_2.0.1   gmodels_2.16.2    purrr_0.2.3      
## [25] xtable_1.8-2      scales_0.5.0.9000 gdata_2.18.0     
## [28] jpeg_0.1-8        ggmap_2.7         tibble_1.3.4     
## [31] withr_2.0.0       lazyeval_0.2.0    proto_1.0.0      
## [34] mime_0.5          deldir_0.1-14     memoise_1.1.0    
## [37] evaluate_0.10.1   nlme_3.1-131      MASS_7.3-47      
## [40] tools_3.3.3       data.table_1.10.4 geosphere_1.5-5  
## [43] RgoogleMaps_1.4.1 stringr_1.2.0     munsell_0.4.3    
## [46] rlang_0.1.2       rjson_0.2.15      htmlwidgets_0.9  
## [49] igraph_1.1.2      crosstalk_1.0.0   bitops_1.0-6     
## [52] base64enc_0.1-3   labeling_0.3      rmarkdown_1.6    
## [55] boot_1.3-19       gtable_0.2.0      reshape2_1.4.2   
## [58] R6_2.2.2          knitr_1.17        bindr_0.1        
## [61] rprojroot_1.2     stringi_1.1.5     Rcpp_0.12.12     
## [64] mapproj_1.2-5     png_0.1-7         coda_0.19-1